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ABSTRACT 

Recent developments in the application of spectral methods to two- 
dimensional compressible flows are reviewed. A brief introduction to spectral 
methods — their history and especially their implementation — is provided. 
The stress is on those techniques relevant to transonic flow computation. The 
spectral multigrid iterative methods are discussed with apfllcatlon to the 
transonic full potential equation. Discontinuous solutions of the Euler 
equations are considered. The key element is the shock fitting technique 
which la briefly explained. 
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1. IirrRODOCTION 

Spectral methods have their roots in approximation 
theory* They are based on representations of the solution to 
a problem by a finite series of global (and preferably ortho- 
gonal) functions* The expansion coefficients are usually re- 
ferred to as spectra, and hence this technique Is called the 
spectral method. All the derivatives of the solution are ap- 
proximated by the corresponding derivatives of the finite 
series expansion. Under the right circumstances such hlgh- 
order approximations can produce extremely accurate numerical 
solutions. There are three versions of spectral methods: 
spectral Galerkln, spectral tau, and spectral collocation. An 
extended discussion of each of these versions Is given in tl]* 

The first serious application of spectral methods to 
fluid dynamics used the Galerkln approach: the solution was 

expanded In a series of functions satisfying the boundary con- 
ditions and the calculation was performed entirely In terras of 
the expansion coefficients. Already in 1954 Sllberraan [2] 
used them for meteorological modelling* The numerous Investi- 
gations which followed Sllberraan^ s work established the feasi- 
bility of the spectral method for low resolution calculations. 
In particular, Ellsaesser [3] showed that for the simple 
balanced baro tropic model the efficiency of a low resolution 
spectral method was competitive with the then available low 
resolution finite difference methods. However, the cost of 
these early spectral Galerkln methods soared to a prohibitive 
level as the number of expansion functions (and hence the 
resolution) Increased. The high cost was due to the straight- 
forward manner In which the convolution sums arising from the 
nonlinear terras were evaluated. 

The breakthrough came when Orszag [4], [5] (and also 

Ellasen, et al* [6]) proposed a transform method for handling 
the conv">lutlon sums ^ This change so Improved the efficiency 




of spectral methods that they became practical for hlgh- 
resolutlon calculations, even in three dimensions. Indeed, 
the accuracy of spectral Galerkin methods is so great, 
especially in terms of their extremely favorable phase errors, 
that they are now routinely used in numerical weather predic- 
tion. They have also been profitably applied to the simula- 
tion of homogeneous, isotropic turbulence [7). 

The spectral tau method was devised by Lanczos [8]. Its 
principal difference from the spectral Galerkin method lies In 
the treatment of boundary conditions. Lanczos used Chebyshev 
polynomials as the expansion functions for solving linear 
ordinary differential equations with rational coefficients. 
Orszag [9], [10] has applied this metnod to certain fluid 
dynamics problems. 

For many problems, especially nonlinear ones, the 
spectral collocation method is the easiest to Implement as 
well as the most efficient. The earliest investigations of 
this method are those of Krelss and Oliger [11], (who called 
it the Fourier method) and Orszag [6], [12] (who termr it 

pseudospectral) . Thus far this has been the only type of spec- 
tral method yet applied to transonic problems. The present 
discussion will be confined to spectral collocation methods, 
with all future references to spectral methods implicitly re- 
ferring to this specific type. 

Spectral calculations of compressible flows have only 
been performed in the last few years. The initial invest iga- 
tious were for one-dimensional flows. These were carried out 
by Zang and Hussalnl [13], Gottlieb, et al. [14], and Taylor, 
et al. [15] . Some promising two-dimensional transonic results 
have been obtained recently for the full potential equation by 
Streett [16] and Streett, et al [17] and for the Euler equa- 
tions by Salas, et al. [18] and Hussaini, et al. [19]. This 
article will describe the details of the spectral methods em- 
ployed in these two-dimensional calculations and will present 
some representative results. Since spectral methods are a 
novel approach to transonic flow computations, a basic intro- 
duction to their properties and implementation will be pre- 
sented first. 


2. FUNDAMENTALS OF SPECTRAL METHODS 

For problems with periodic boundary conditions spectral 
\ methods based upon Fourier series are the obvious choice. If 

the boundary conditions are Dirlchlet or Neumann, then Cheby- 
shev polynomials are usually employed. For many problems an 
appropriate spectral method can produce a rapidly convergent 
approximation. The particular choices mentioned above have 
the added advantage of efficient Implementation via the Fast 
Fourier Transform. This section will furnish the details 
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necessary for implementing these specttal methods* To begin 
with, however, the convergence properties of spectral methods 
will be illustrated with two eleme tary examples. The first 
of these has some of the characteristics of the steady poten- 
tial calculations described in Section 4, whereas the second 
is relevant to the time-dependent Ei^ler solutions of Section 
5. 


2.1. A One-dimensional Elliptic Problem 

A major reason for the appeal of spectral methods is 
their potential for rapidly convergent approximations. The 
Fourier series expansion over [0,2tt] of the function 


u(x) 


3 

5-4 cos X 


( 1 ) 


provides a simple Illustration. It is 


u(x) =» I u e , 


where 





( 2 ) 


(3) 


The approximation obtained by truncating the Fourier series, 

N/2.-1 


u„(x) « ^ 

^ k— N/2+1 '' 


\ ~ ikx 

L u. e , 


(4) 


satisfies 


|Uj^(x) - u(x) I < 4e 




(5) 


It converges exponentially fast, i.e., 

N^|u„(x) - u(x)| 0 as N -► 00 (6) 

N 


for all positive Integers p. This property of exponential 
convergence is exhibited by the truncated Fourier series of 
\ any periodic function which is infinitely differentiable* The 

rapid decay of the Fourier series of such a function follows 
from repeated Integrations-by-parts: Let u(x) be periodic 

and infinitely differentiable on [0,2tt], Its Fourier coef- 
‘ ; ficlents are given by 
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1 


I , , -Ikx, 

} u(x)e dx. 


(7) 


0 


A single Integration-by-parts yields 

2tt 
0 


~ _ 1 
"k (-ik) 


-Ikx 


|2tt 2Tf 

u(x) 1 + TTlTr / u'(x)e ’^dx. 

In Q 


( 8 ) 


The boundary term vanishes because of the periodicity require- 
ment. The remaining integral is o(l) as k -► ® by the 
Riemann-Lebesque Lemma because of the differentiability condi- 
tion. Thus, a single integration lets us conclude that 

" o(k"^). Clearly, p Integrations produce - o(k'^P). 
Note that both the periodicity and differentiability condi- 
tions are necessary for this argument. If either fails, then 
the rate of convergence Is algebraic. 

A Fourier spectral method for a differential equation 
makes use of some finite Fourier series representation of the 
solution. This series will be related to, but different from, 
the truncated Fourier series of Eq. (4). The details are fur- 
nished in Section 2.3. Consider the differential equation 



(9) 


on [0,2tt] with periodic boundary conditions. Suppose that 

f is chosen so that the exact solution is given by Eq (1). 

Applied to this problem the spectral method yields the results 
shown in the third column of Table I. The second column gives 
the results of the truncated Fourier series and the last col- 
umn reports the results for a second-order finite difference 
method. Note that the bound given by Eq. (5) for the trunca- 
ted Fourier series is sharp — the entries In the second col- 
umn agree precisely with this bound until N Is so large that 
round-off error predominates. (These calculations were per- 
formed on a CDC Cyber-175, which has 14 significant digits.) 
The correct N * 128 and N = 256 entries are 2. 17(-19) and 
1.18(-38), rerpectlvely. The exponential convergence of the 
truncated Fourier series is evident. The spectral method per- 
forms nearly as well. (In fact. Its maximum error for N " 32 
and N ■ 64 is considerably lower, but this Is due to a can- 
cellation of error. In an RMS sense these two respective 
truncated series and Fourier spectral errors are within sever- 
al percent.) Only at very low resolution is the spectral 

method substantially worse than the truncated series. But, of 
course, the major point of this example is the decided superi- 
ority of the spectral method over the second-order finite dif- 
ference approximation. 


f 
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Table !• HexlauB Error for a 1*^ Periodic Problem 


N 

Truncated 

Series 

Fourier 

Spectral 

Finite 

Difference 

4 

1.00 

(0) 

4.42 

(0; 

6.28 

(0) 

8 

2.50 

(-1) 

1.28 

(0) 

1.98 

(0) 

16 

1.56 

(-2) 

3.15 

(-2) 

2.02 

(-1) 

32 

6.10 

(-5) 

3.46 

(-5) 

3.61 

(-2) 

64 

9.31 

(-10) 

4.89 

(-10) 

8.65 

(-3) 

128 

5.68 

(-14) 

7.11 

(-14) 

2.14 

(-3) 

256 

5.68 

(-14) 

7.11 

(-14) 

5.34 

(-4) 


R. f 


2, 2. A One-dimensional Hyperbolic Problem 


Spectral methods for time-dependent problems can also ex- 
hibit exponential convergence. Indeed, spectral methods have 
thus far made a greater Impact on evolution equations than on 
steady-state ones* A simple example is provided by the wave 
equation 


3u . 

3t 9x 


0 


( 10 ) 


on the interval [-1,1] with initial condition u(x,0) and 

boundary condition u(-l,t)* Since this is not a periodic 
problem, a spectral method based upon Fourier series in x 
would exhibit extremely slow convergence: the Fourier coeffi- 
cients decay only as fast as 0(k"^) in the general case* 
(The integration-by-parts argument given earlier falls, even 
if the solution is infinitely differentiable because the 
boundary term in Eq. (8) is finite.) Hov^ever, rapid conver- 
gence as well as efficient algorithms can be attained for 
spectral methods based upon Chebyshev polynomials. These are 
defined on [-1,1] by 

T (x) “ cos (n cos ^x). (11) 


The f unc t ion 

u(x,t) * sin otTT(x-t) (12) 


is one solution to Eq. (10). It has the Chebyshev expansion 

00 

u(x,t) - I Z(t) T (x), (13) 

n-0 " 


where 
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u (t) - — sin - aiTt) J (ot) 
n c 4 n 

n 


n * 0 
n > 1 


and Bessel function of order n. The asympto- 

tic properties of the Bessel functions imply that 


n*^u (t) ^ C 


as n ^ 


for all positive integers p. Note that this result holds 
whether or not ot is an integer. In contrast, the Fourier 
coefficients of u(x,t) are 


sin (oH-k) 1 -laTTt sin (a*k) 
“ 2TT ® OH-k " 2TT ® a-k 


For non-integer ot these decay extremely slowly. 


The change of variables 


X « cos 


the definition 


v(0,t) “ u(co8 6,t) , 


and Eq. (11) reduce Eq. (13) to 


v(9,t) 


u (t) cos n^^ 
n 


Thus, the (Jiebyshev coefficients of u(x,t) are precisely the 
Fourier coefficients of v(0,t). This new function is auto- 
matically periodic# If u(x,t) is infinitely differentiable 
(in x) , then v(0,t) will be infinitely differentiable (in 6). 
Hence, straightforward integral lon-by-parts arguments lead to 
the conclusion that the Chebyshev coefficients of an infinite- 
ly differentiable function will decay exponentially fast. 
Note that this holds regardless of the boundary conditions. 

The effectiveness of this type of approximation is demon- 
strated by the results in Table II. Shown there are the 
errors at t - I for several approximations to Eq. (10) with 
initial and boundary conditions based on Eq. (12) for a - 2.5. 
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(The explicit time-step used In these calculations has been 
taken to be so small that time-disc retlzatlon errors, but not 
round-off errors, are negligible.) The results of a spectral 
ra-ithod based on Fourier series have been included to emphasize 
the Importance of the proper choice of expansion functions. 


Table II* Maxlnm Error for a 1-D HTperbollc Problem 


N 

Truncated 

Series 

Chebyshev 

Spectral 

Fourier 

Spectral 

Finite 

Difference 

4 

1.24 (0) 

1.49 (0) 

1.85 (0) 

1.64 (0) 

8 

1.25 (-1) 

6.92 (-1) 

1.92 (0) 

1.73 (0) 

16 

7.03 (-6) 

1.50 (-4) 

2.27 (0) 

1.23 (0) 

32 

1.62 (-13) 

3.45 (-11) 

2.28 (0) 

3.34 (-1) 

64 

1.79 (-13) 

9.55 (-11) 

2.27 (0) 

8.44 (-2) 


2.3. Implementation of Fourier Spectral Methods 

The key to the implementation of Fourier spectral methods 
Is the Discrete Fourier Transform. Let a function u(x) be 
represented by its values u. « u(x^) at the collocation 
points 

Xj=-^ j - 0,1, (21) 


The discrete Fourier coefficients of Uj are 

“k*N ^ ^ 

j-0 

The Inverse relationship Is 
N/2 - 1 ^ 

Uj * u^e J J = 0,1, •••,«-!, 

k— N/2 


and the orthogonall" / relation Is 


1 

N 


N-1 

I 


j-0 


Ikx j 
e J 



» 


( 22 ) 


(23) 


(24) 


k,l 


where 


1 

0 


k - i, £±N, 1±2N 
otherwise 


(25) 
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Note the difference between Eq. (24) and the orthogonality re- 
lation fo\- continuous Fourier series* In the latter case the 
sum is i‘v 2 pi aced by an Integral and the usual Kronecker delta 
function appears on the right-hand side* 

The counterpart of the periodicity condition 
u(x+2^) * u(x) is 


u 


J+N 


(26) 


Unlike the continuous case, discrete Fourier series also have 
a periodicity requirement on the coefficients: 


^k+N * 


(27) 


This is an immediate consequence of Kq* (22). The connection 
between the discrete and continuous Fourier coefficients fol- 
lows from Eqs. (2), (22), and (24) and is 

00 

^ < 28 ) 

£a-00 

All but the £ » 0 contribution to the sum are the aliases of 
u|^, i.e., Fourier components which are indistinguishable from 
u]^ on the discrete grid. These are a source of error in 

spectral methods in addition to the error that arises from the 
truncation of the exact Fourier series. It is just such 
aliasing terms that account for the differences between the 
entries in columns 2 and 3 of Table T. 


Here then are the details of the spectral method used for 
the Table I results: 


1) Given fj - ^(^j) J * find the dis- 
crete Fourier coefficients for k » -N/2, 

- N/2+1, • • • ,N/2-l . (Use Eq. (22) with f in place 
of u. ) 


2 ) 



|k| < N/2 
k - -N/2 


(29) 


(This is the solution of Eq. (9) in terms of the 
Fourier coefficients.) 


3) Compute uj Itself for j - 0,1,***,N-1 from its 
discrete Fourier coefficients. (Use Eq. (23).) 
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The Kast Fourier Transform (FFT) is an efficient and 
widely-avallable algorithm (even in assembly language) for ac- 
complishing seeps (1) and (3). A broad survey of FFT's is 
provided in [20] • The only slight complication is that stan- 
dard versions of the FFT (such as the IMSL subroutine FFT2) 
produce the coefficients for k * The periodic- 
ity relation given by Eq . (27) enables the desired coeffi- 

cients to be extracted readily from the FFT output. The pe^^ " - 
odlcity relation also helps in using the FFT to perform 
such as those in Eq. (23). 

The ease with which a direct solution is attainable for 
this spectral discretization of Eq. (9) is exceptional. For 
example, no efficient direct solution scheme exists for the 
spectral solution of 

1^ (a(x,y) -^) + Ij;: (b(x,y) - f. (30) 

One must, in general, resort to iterative schemes for its so- 
lution. An essential element of these schemes is the explicit 
evaluation of terms such as those appearing in Eq. (30). Con- 
sider just the first term. Given u(x,y) at the collocation 
points, this terra is evaluated by 

1) computing 3u/3x by Fourier collocation, 

2) multiplying by a(x,y), and 

3) computing 9/9x(a(x,y) 8u/9x] by Fourier collocation* 

The differentiation occurring in step 1 Involves (with the y- 
dependence suppressed): 

(i) using the FFT to compute the discrete Fourier coeffi- 
cients u]^ for k » -N/2,-N/2+l, • •• ,N/2-l, 

( 11 ) 

( 111 ) 


The choice of calls for some explanation. Since 

u(x) ^ Is real, so too will be But only the real part 

of v«^/2 roakes any effective contribution to 9u/9x. If a 
derivative of higher-order Is desired, then steps (1) to (111) 
still apply but with the appropriate power of Ik appearing 
as the multiplier of In Eq. (31). 

The FFT enables the le -ha^u side of Eq. (30) to be 
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evaluated at all the N x N collocation points In 0(N^^n N) 
operations* For large N this compares unfavorably with the 
finite difference cost of O(N^) for an N ^ N grid. However, 
when judged by the cost of equally accurate solutions the 
spectral method Is likely to be cheaper since it needs a 
coarser grid (sec Table I)* 

2*4. Itnplementatl> " n of Chebyshev Spectral Methods 

A Chebyshev spectral method makes use of c finite 
Chebyshev series such as 


N 

u^(x) = I u^T^(x). (32) 

n=»0 

The standard collocation points are 


= cos 


j « (33) 


Thus , 


where u. 
tlon Is 


N 

Uj = ^ 

n^O 

is the approximation to u(x^). The Inverse rela- 


where 


# } "i"' “,i ’ 

”^n j“0 



j =■ 0 or N 
1 j ^ N-1 


(35) 


(36) 


These last two sums may be evaluated by the FFT* Tht standard 
FFT, however, does complex arithmetic and Ignores the symmetry 
present In a cosine transforni. The second appendix of [6] 
describes how to make more efficient use of the FFT for 
evaluating the sums In Eqs. (34) and (35). One may also use a 
Fast Cosine Transform. A Fortran listing of one version Is 
given in [21]. 


The Chebv.aev collocation points are the extreme points 
of T^(x). Note that they are not evenly distributed In x, 

but rather are clustered near the endpoints* The smalleot 
mesh size scales as 1/N^* While this distribution contri- 
butes to the quality of the Chebyshev approximation and per- 
mits the use of the FFT in evaluating the series. It also 




u 


C;v. . . 

OP f'Oun g 


r 


TY 


places a severe time-step limitation on e.-.pllclt rjethods for 
evolution equations. 

A Chebyshev spectral method for Eq. (10) combines some 
explicit tlme-dlscretlzatlon with an approximation to the -spa- 
tial derivative which Is based upon analytical differentiation 
of the Chebyshev series for u. Consider first the Infinite 
series, for which 

00 

u(x) => I u T (x), (37) 

n n n 


with the tliTie dependence of u suppressed. Write the expan- 
sion of the derivative as 


I 


u'(x) 




(x) . 


( 38 ) 


The goal Is to relate the coefficients In Eqs. (37) and 
(38). The starting point Is the recursion relation 


iH-l n-1 


_2 

c 

n 


T (x), 
n 


(39) 


which follows from Eq. (11). Inserted into Eq. (38) this pro- 
duces 


« T' , (x) T' , (x) 

^ \ 1 / V (1 ) r n+1 n-1 

u (x) I c u [ r; T“ 

n»0 ” ^ 


^ n-1 n-1 


:(i) 


n»l 

But from Eq. (37) 


2n n 


T'(x) - )' -?^T:(x). 


n*l 


2n n 


u (x) “ / u T (x) . 

' - n n 
n=^l 


Therefore, 


2n u . c 

n n-1 n-1 n+1 


n > 1. 


(40) 


( 41 ) 


( 42 ) 


The evaluation of the 9u/3x term In the Chebyshev spec- 
tral method for Eq. (10) consists of: 

1) Given Uj ■ u(xj) for j ■ 0,^1,*»*,N find the dis- 
crete Chebyshev coefficients u^ for n - 0, !,•••, N. 
(Use Eq. (35).) 


n 


r 







3) Compute 9u/^x(Xj) Itself for j « 0, !,•••, N from 
its discrete Chebyshev coefficients. (Use the analog 
of Eq . (34) .) 

Higher-order derivatives can be calculated by repeating step 
(2) as often as needed* When the FFT is used for steps (1) 
and (3) the cost of a derivative evaluation is 0(N £n N). 

Eor the Eq. (10) calculation, the derivative 9u/9x is 
not needed at the inflow boundary (x « -1) since the boundary 
condition rather than the PT)E is used to update u^. Note 

that there is no need for a special formula at the outflow 

boundary (x - +!)• Although the PDE is used to update Uq, 

the value of 3u/9x at x ■ 1 is automatically available 

from the Chebyshev spectral calculation outlined above. In 
concrast the second-order finite difference calculation used 
for Table II employed the special formula of first-order 

extrapolation at the outflow boundary. 

As a general rule the correct numerical boundary condi- 
tions for a spectral method are the same as the correct analy- 
tical boundary conditions. The global nature of the approxi- 
mation avoids the need for special differentiation formulae at 
boundaries. At the same time spectral methods are quite un- 
forgiving of Incorrect boundary ccndi.lons. The inherent dis- 
sipation of these methods is so low that boundary errors 
quickly contaminate the entire solution. In many fluid dyna- 
mical applications the computational region must be terminated 
at some finite, artificial boundary. The difficulty at 
"artificial” boundaries is that analytically correct, fully 
nonlinear boundary conditions for systems are seldom known. 
One example of a workable artificial boundary condition for 
the Euler equations is given Ir Section 5.4. 


I 


3. SPECTRAL MDLTI6RID ITERATIVE METHODS 


Direct solutions of implicit spectral equations 
wh -her arising from elliptic problems or from implicit time- 
discretizations of evolution problems — are rarely feasible. 
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Except in special cases the matrices representing spectral 
apf roximatlons are full. Iterative methods are a practical 
alternative because the requisite matrix-vector products can 
be evaluated via the FFT. 

3.1. An Elementary Example 

An attractive iterative scheme for spectral equations 
utilizes multigrid concepts. The basic description of spec- 
tral multigrid methods for linear, elliptic equations is given 
in [ 22 ] and [23]. Additional considerations for the non- 
linear, potential flow application are given in [17]. A brief 
summary of these concepts is given here since they play a 
major role in the spectral transonic potential flow calcula- 
tions discussed in the following section. 

The fundametitals of spectral multigrid are perhaps 
easiest to grasp for the simple model problem 



on [0,2 tt] with periodic boundary conditions. The Fourier 
approximation to the left-hand side of Eq. (45) at the collo- 
cation points is 


P= 


N/2-1 

I 


-N/2+1 




Ipx. 


(46) 


The spectral approximation to Eq. (45) may be expressed as 


LU - F, (47) 

where 

U - (uq,u^,*“»Uj^__j), (48) 

F - (fQ,fp--*,fjj_,j), (49) 

and L represents the Fourier spectral approximation to 

- d^/dx^. 

A Richardson's iterative scheme for solving Bq. (47) Is 
V V + w(F - LV), (50) 



■i 
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where <*) Is a relaxation parameter, on the right side of the 
replacement s 3 nnbol (^) V represents the current approximation 
to U, and on the left It represents the updated approximation. 
The eigenfunctions of L are 

£. , . 2irijnAi .cix 

Cj(p) » e , (51) 

with the corresponding eigenvalues 

A(p) - p^ (52) 

where j « and p - N/2+1, • • • ,N/2-le The 

Index p has a natural interpretation as the frequency of the 
eigenfunction* 

The error at any stage of the iterative process is V - U; 
it can be resolved into an expansion in the eigenvectors of 
L* Each iteration reduces the p'th error component to v(A ) 
times its previous value, where ^ 


v(X) =1 1 — ujA, 


(53) 


The optimal choice of results from minimizing |v(X)) for 

^ ' t^min^ ^max^» \±n " ^ ^max “ 

need not worry about the p ® 0 eigenfunction since It corre- 
sponds to the mean level of the solution, which is at one's 
disposal for this problem*) The optimal relaxation parameter 
for this single-grid procedure is 


2 

. 

SG X + X ^ 
max min 

It produces the spectral radius 


(54) 


X ^ X ^ 

P ^ max min 

“ X + X ^ • 

max min 

Unfortunately, - 1 - 8/N^ , which implies that 
iterations are required to achieve convergence* 


(55) 

0(n2) 


This slow convergence is the outcome of balancing the 
damping of the lowest-frequency eigenfunction with that of the 
hlghest-frequency one In the mlnlmax problem described after 
Eq. (53). The multigrid approach takes advantage of the fact 
that the low-frequency modes (|p| < N/4) can be represented 
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just as well on coaiset* ^rids. It settles for balancing the 
middle-frequency eigenfunction (|p| = N/4) with the highest- 
frequency one (|p| = N/2), and hence damps effectively only 
those modes which cannot be resol /ed on coarser grids. In 
Eqs. (54) and (55), ^j^ln replaced with • 

The optimal relaxation parameter In this context Is 


0) 


MG 


X + X _ 

max mio 


The multigrid smoothing factor 


^MG 


X - X , - 
max mid 

max mid 


(56) 


(57) 


measures the damping rate of the high-frequency modes. In 
this example Mwq ** 0.60, independent of N. The price of 
this effective damping of the high-frequency errors Is that 
the low-frequency errors are hardly damped at all* However, 
on a grid with N/2 collocation points, the modes for 
Ipl e [n/ 8, N/4] are now the high-frequency ones. They get 
damped on this grid. Still coarser grids can be used until 
relaxations are so cheap that one can afford to damp all the 
remaining modes, or even to solve the discrete equations 
exactly. 


Let us consider just the interplay between two grids. A 
general, nonlinear fine-grid problem can be written 


L^(U^) - F^. 


(58) 


The shift to the coarse grid occurs after the fine-grid 
approximation has been sufficiently smoothed by the 

relaxation process, l.e., after the high-frequency content of 
the error has been sufficiently reduced. The 

related coarse-grid problem Is 


L^(U^) » F^, 


(59) 


where 


- r[f^ - L^(V^)] + l'^(RV^). ( 60 ) 


The restriction operator R Interpolates a function from the 
fine grid to the coarse grid. The coarse-grld operator and 
solution are denoted by L^ and U^, respectively. After an 
adequate approximation to the coarse-grid problem has 
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been obtained, the fine-grid approximation Is corrected via 


yf ♦ + P(V<^ - RV^). (61) 


The prolongation operator P interpolates a function from the 
coarse grid to the fine grid. 

3.2. Interpolation Operators 

The spectral raultigrld interpolation operators for perio- 
dic coordinates amount to trigonometric interpolation: for 

example, given a function on a coarse grid, compute the dis- 
crete Fourier coefficients and then use the resulting discrete 
Fourier series to construct the interpolated function on the 
fine grid. This may be accomplished by performing two FFT^s. 
Interpolation for non-periodic coordinates employs Chebyshev 
series in an analogous fashion. Detallc:d descriptions of the 
interpolation operators are available in [17] and [23]. 

3.3. Coarse-Grid Operator 

A typical term in the class of problems containing 
potential flow is 

[a(u,x) ^]. (62) 

The discrete operator Which represents its fine-grid spectral 
approximation Is 

- 0 A 0, (63) 


where 0 is a spectral first-derivative operator (either 
Fourier or Chebyshev) and A is the diagonal matrix 




(64) 


Many multigrid investigators have advocated choosing the 
coarse-grld operator so that 


- RL^P. 


(65) 


Both the Fourier and the Chebyshev first-derivative operators 
satisfy 


qc • RO^P. 


( 66 ) 


r 



However, Eq. (65) Itself is not satisfied if the coarse-grid 
analog of Eq. (63) is used to define L^, except in the 
trivial case for which a(u,x) is a constant. On the other 
hand, much of the efficiency of the pseudospectral method Is 
lost if Eq. (65) is used to define the coarse-grid operator. 
The most satisfactory compromise seems to be using Eq. (63) 
but with the restricted values of a(u^,Xj) in place of the 
pointwise values. 


4. THE COK?RESSIBLE FOTEIITIAL EQUATION 

The computational results of the past decade have demon- 
strated that fairly accurate predictions for a number of tran- 
sonic flows of practical Interest can be made on the basis of 
the compressible potential equation. This nonlinear equation 
is of mixed elliptic-hyperbolic type, precluding purely ellip- 
tic or purely hyperbolic solution procedures. The numerical 
solution of the potential equation became feasible only after 
the introduction of type-dependent differencing by Hirman and 
Cole [24]. The review by Hall [25] provides an exhaustive 
history of computational approaches to the potential equation. 

Until the recent work of Streett [16] , the discretization 
procedures for the potential equation were invariably based on 
low-order finite difference or finite element methods. 
Streett used a spectral discretization of the full potential 
equation and obtained its solution by a single-grid iterative 
technique* The application of spectral multigrid techniques 
by Streett, et al. [17] produced a dramatic acceleration of 
the iterative scheme. Even in its relatively primitive state 
the spectral multigrid scheme is competitive, and in some 
cases unequivocally more efficient, than standard finite dif- 
ference schemes. 

4. 1. The Reduced Potential Problem 


Streett solved the two-dimensional full potential equa- 
tion (applying boundary conditions at the actual airfoil sur- 
face) • In this work a numerical conformal mapping (also gen- 
erated by Fourier techniques) was used to transform the air- 
foil onto the unit circle. Moreover, the calculations were 
actually performed in terras of the reduced potential G, which 
is defined by 


^ “ (R + r) cos .0 - E tan ^ tan o], 


(67) 


where ^ is the potential, R and 0 are the computational 
polar coordinates, E is the circulation and is the Mach 

number at infinity. The reduced potential is periodic in 0 
and it satisfies 
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iG' . 3_f£ in') 
9R'' 30'‘R ae-* 


0 , 


along with 



G + 0 


at R ” I, 


as R + 


( 68 ) 


(69) 

(70) 


and the Kutta condition* The density Is given by the isentro- 
pic relation ^ 

+ 4 - 1 ) 1 '^; <”) 

the ratio of specific heats Is denoted by Y, the velocity 
components In the physical (r,9) plane are 


i ii 

**r “ H 


(72) 


“ RH 90 ’ 


(73) 


4 0 

and the Jacobian between the complex physical plane (z ■ re^ ) 
and the complex computational plane (o Re*’*^) Is 


H 



(74) 


4.2. Discretization 


The spectral method employs a Fourier series representa- 
tion in 0. Constant grid spacing in 0 corresponds to a 
convenient dense spacing In the physical plane at the leading 
and trailing edges. The domain in R (with a large, but 
finite outer cutoff) is mapped onto the standard Chebyshev do- 
main [-1,1] by an analytical stretching transformation that 
clusters the collocation points near the airfoil surface. The 
stretching is so severe that the ratio of the largest-to- 
smallest radial Intervals is typically greater than 1000. 

The most expedient technique for dealing with the mixed 
elliptic-hyperbolic nature of the transonic problem la to use 
the artificial density approach of Hafez, et al. [26]. The 
original artificial density Is 


p - p - y6p 


(75) 


t 
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M 


(76) 


where M Is the local Mach number and Is an upwind 

first-order (undivided) difference. The spectral calculations 
employed a higher-order artificial density formula. The 
spectral method also required a weak filtering technique to 
deal with some high-frequency oscillations generated by the 
shock. Details are available in 116]. 


4.3. Spectral Multigrid Solution Scheme 

Let the spectral discretization of Eqs. (68) - (70) be 
denoted by 

M(U) = 0, (77) 

and define 

J(U)-|^(U). (78) 


f 


A suitable relaxation scheme for the spectral multigrid solu- 
tion of transonic potential flow is based upon approximate 
factorization techniques similar to those used in finite dif- 
ference discretizations [27]* The Jacobian J(U) is split 
into the sum of two operators ^^^ch invol- 

ving derivatives in only the one coordinate airection indica- 
ted by the subscript. The most straightforward spectral ap- 
proximate factorization scheme is 


[al - J (V)] [al - J (V)]AV « o)aM(V) , (79) 

X y 

where V is the last estimate of U, and V + AV is the next 
estimate. This is commonly referred to as AFl for the tran- 
sonic problem. For second-order spatial discretizations the 
term [al - J^(V)] leads to a set of trldlagonal systems, one 
for each value of y. The second left-hand side factor pro- 
duces another set of tridiagonal systems. For spectral dis- 
cretizations, however, these systems are full; hence, Eq. (79) 
is still relatively expensive to Invert. The spectral factor- 
ization scheme makes the additional approximation of replacing 
J^ and J with their second-order finite difference 
analogs, denoted by and Hy, respectively: 

[al - H (V)] [al - H (V)]AV - u)aM(V). 

X y 


if 


(80) 
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For purely finite difference approximations some analy- 
tical results are available for selecting optimal values for 
the parameters ot and o) [28]. No similar results are yet 
available for the spectral approximation. By analogy with the 
finite difference case was chosen to be of order unity and 
a senuence of ex's was selected in a range 

a - . (8‘) 

where K denotes the number of distinct ot^s. The choices 
of and were based in part on estimates ot the eigen- 

value range of the discrete operators and in (much greater) 
part by trial and error- Fortunately, the AFl scheme is not 
very sensitive to these parameters. This basic iterative 
scheme may be employed in either a single-grid or a multigrid 
context. In the latter case the parameters and 

should be chosen separately for each grid to optimize the 
high-frequency damping. 

The spectral raultigrid solutions of Streett, et al., used 
three different fine grids (with the coarser levels in 
parentheses): 16 x 32 (12 x 16 and 8 x 8), 16 x 48 (14 x 2 ?.^ 

12 X 16 and 3x8) and 18 x 64 (16 x 48 , 14 x 32 , 12 x 16 and 
8x8). Note that in passing to a coarser level the grid is 
typically reduced by less than a factor of 2 in each 
coordinate direction. This choice leads to a significant 
Improvement over the standard gridding for the spectral 
potential flow problem, especially in the supercritical regime 
where the solution has large high-frequency content. 

All the spectral multigrid results were obtained with the 
same fixed schedule: start on the finest grid, work down to 

the coarsest grid and th^-n back up to the finest grid; on the 
way down there is \ sweep though the (three) parameter 
sequence and on the way up there are 2 sweeps. 

4 . 4 . Airfoil Exaraples 

The flow past an NACA 0012 airfoil at 4® angle of attack 
and a freestrearo Mach number of 0.5 is a challenging subcrl- 
tlcal case. rhe airfoil produces a fairly large lift coeffi- 
cient at these conditions and the surface pressure distribu- 
tion shows a sharp suction peak near the leading edge. Since 
the local Mach number in this peak is nearly 1, compressibili- 
ty effects are substanclal. 

Nevertheless, the spectral solution on a relatively 
coarse grid captures all the essential details of the flew. 
The surface pressure coefficient from the spectral code MGAFSP 
[17] using 16 points in the radial (R) direction, and 32 
points in the azimuthal (0) direction is displayed in Fig. 1. 
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Figure 1. Spectral (triangles) and finite difference (solid 
line) surface pressures for a suberic leal floir. 

The symbols denote the solution at the collocation points* 
For comparison, the result from the finite difference, multi- 
grid, approximate factorization code FL036 [29] is shown as a 
solid line. The grid used in the benchmark finite difference 
calculation is so fine (64 384 points) that the truncation 

error is well below plotting accuracy. The FL036 and MGAFSP 
results are identical to plotting accuracy. The spectral 
computation on this mesh yields a lift coefficient with 
truncation error less than 10***^. Spectral solutions on a 
16 X 32 grid are thus of more than adequate resolution and 
accuracy for subcrltical flows. 

In Figure 2 are shown convergence histories from FL036, 
MGAFSP, and the finite difference, approximate factorization, 
single-grid code TAIR [27]. Meshes which yield approximately 
equivalent accuracy were chosen. The surface pressure results 
are the same , to plotting accuracy, the lift coefficient 1s 
converged in the third decimal place, and the predicted drag 
coefficient is less than .001. (Actually, the spectral result 
Is an order of magnitude more accurate than these limits, uut 
the TAIR result barely meets them.) 

A lifting sup*-:rcrltlcal case is provided by the NACA 0012 
airfoil at - 0.75 and a ■ 2®. This yields a section 

lift coefficient of nearly O.6. A shock appears only on the 
upper surface for these conditions and Is rather strong for a 
potential calculation; the normal Mach number ahead of the 


r 
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Figure 2. Nulmm residual versus nachlne 
time for a subcrltlcal flov* 
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Tlgttre 3* Surface pressures for a supercritical flow: 

a) MfiAFSP at collocation points; b) MCAFSP Inter- 
polated onto finer grid; c) TAIR, and d) PL036. 
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shock is about 1.36* Lifting supercritical cases are 
especially difficult for spectral methods since the solution 
will always have significant content in the entire frequency 
spectrum: tVie shock populates the highest frequencies of the 

grid and the lift is predominantly on the scale of the entire 
domain. An iterative scheme therefore must be able to damp 
error components across the spectrum. 

Surface pressure distributions from MOAFSP, TAIR, and 
FL036 are shown in Figure 3. The respective computational 
grids are 18 x 64, 30 x 149 ^ and 32 x 192 , The latter two 
are the default grids for the production finite aifference 
codes. Spectral results obtained by trigonometrically inter- 
polating the 18 X 64 grid results onto a much finer grid are 
included alongside the results at the collocation points. 
This reveals the wealth of detail that is provided by the 
rather coarse spectral grid. The shock predicted by TAIR is 
far more rounded and smeared than that of FL036, reflecting 
the coarser mesh and larger artificial viscosity used in the 
former. The TAIR result shown is also only correct to one 
decimal place In lift as compared with a flner-grld result. 
Convergence histories for these three cases are shown in 
Figure 4 along with the results for MCAFSP on a coarser grid 

(16 X 48). 



Figure 4. NaxlMum residual versus aachloe 
tlae for a supercritical floir. 

5. THE EULER EQUATIOHS 

The Euler equations undoubtedly provide more information 
than the potential flow approximation. The numerical diffi- 
culties in solving the filler equations are well known. The 
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problems tend to be even more severe wit-i spectral than finite 
difference methods. Explicit time-stepping schemes are 
especially costly because the Chebyshev collocation points 
have a very small spacing near the boundary. The numerical 
boundary conditions, particularly for artificial boundaries, 
must be sophisticated because spectral methods are extremely 
sensitive to Improper boundary treatments. The oscillations 
arising In shock-capturing methods are quite troublesome 
because the global nature of spectral approximations spreads 
the oscillations over the entire domain. 


Spectral methods for compressible flows are still so 
novel that most of these difficulties remain to be surmounted. 
However, at least the shock-induced oscillations can be 
avoided by resorting to shock-fitting techniques. Here the 
shock front is a computational boundary whose shape and motion 
are generated during the calculation. Since the flow within 
the computational domain is smooth, there la reason to expect 
a shock-fitted solution to be highly accurate. Shock-fitted 
spectral solutions to the r.uier equations were first presented 
by Salas, et al. [18]. Results for related problems were sub- 
sequently given by Zang, et al. [30]. Additional examples and 
more numerical details are contained In [19]. The essential 
features of these Investigations follow. 

5« 1* The Shock Interaction Problem 
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Figure 5. Ifodel problea In the pJtyslcal dowln e short 
tlse after the start of the calculattoa. 
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A model problem which has been used to study the 
interaction of a shock wave with a vortex [31] or with 

idealized turbulence [32] is Illustrated in Figure 5, At 
time t = 0 an infinite, normal shock at x ** 0 separates a 
rapidly moving, uniform fluid on the left from the fluid on 

the right which is in a quiescent state except for some 
specified fluctuation. The initial conditions are chosen so 
that in the absence of any fluctuation the shock moves 
uniformly in the positive x-dlrectlon with a Mach number 
(relative to the fluid on the right) denoted by M^. In the 
presence of fluctuations the shock front will develop ripples. 
The shape o^ the shock is described by the function Xg(y,t'. 

The numerical calculations are used to determine the state of 

the fluid in the region between the shock front and some 
suitable left boundary slso to determin the 

motion and shape of the shock front itself. 


Tlie physical 
is given by 


domain in which the fluid motion Is computed 

X (t) < X < X (y,t) 

L s 


( 


—00 < y < «> 


( 82 ) 


t > 0. 


The change of variables 

X - Xj^(t) 

^ “ Xg(y.t) - Xj^(t) 

Y - I [l + tanh(ay) ] (83) 

T - t, 


produces the computational domain 


0 < X ^ 1 

0 < Y < 1 (84) 

T > 0. 

The stretching parameter a is typical of order 1. 

The fluid motion is modeled by the two-dimensional Euler 
equations. In terras )f the computational coordinates these 
are 




r 
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whore 


Q = tP,u,v,S] , 


B = 



and 



U 

0 

0 


V 

0 

0 


0 

u 

0 


0 

V 

0 


0 

0 

0 

U ^ 

0 "■ 
0 

0 

V 


The contravarlant velocity components are given by 


( 85 ) 


( 86 ) 


(37) 


( 88 ) 


U - X 

+ uX 

+ vX 

t 

X 

y 

V = Y 

+ uY 

+ vY . 

t 

X 

y 


(89) 


A subscript denotes partial differentiation with respect to 
the indicated variable. P, a, and S are the natural 
logarithm of pressure, the sound speed, and the entropy 
(divided by the specific heat at constant volume), 
respect!' ‘^ly, all normalized by reference conditions at 
downstream infinity; u and v are velocity components in 
the X- and y-directions, both scaled by the characteristic 
velocity defined by the square root of the pressure-density 
ratio at downstream infinity* A value Y = 1*4 has been 
used* 


5.2. Dj scretizatio n 

Let k denote the time level and let At be the time- 
step increment. The time discretization of Eq. (85) is then 
as follows; 


Q 


1 - AtL^jq' 


kl^k 


» 


(90) 
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[Q*" + (1 - AtL)0], (91) 

where the spatial operator L represents an approximation 
to B 9/3x + C 9/9Y. In the spectral method, the solution 
Q is first expanded as a double Chebyshev series, 

M N 

Q(X,Y,T) = I I Q (T)T (C)T (n), (92) 

p.O q-0 P9 P 9 

where 

C « 2X - 1 and n * 2Y - 1, (93) 


and Tp and Tq are the Caebyshev polynomials of degrees 
p and q. The derivatives appearing in the spatial operators 
are then evaluated as 

M N ,, ... 

Qy = 2 I I 0^„’"X<^>T^(n), (94) 

^ P -0 q «0 P‘1 P *1 

and 

M M 

Q =2 f I (?)T^(n). (95) 

Y p^o pq P q 

The Chebyshev coefficients of the X-derlvative are denoted by 

They are evaluated by the recursion formulae of Eos* 

(A3) and (44) for each q. The Y-derivatlve is handled in a 
similar fashion. 

Spectral methods for all but constant-coefficient, linear 
problems require some sort of weak filtering for stability. 

For the calculations presented below, the upper third of the 
frequency spectrum of the lolution was filtered every 50 time- 
steps or so. Details are given in [19]. 

5.3. Shock Fitting 

The Ranklne-Hugonlot conditions are used both to 
determine the flow variables (P, u, v, and S) immediately 

upstream of the shock and to determine the shock position. 
Use the subscripts I and 2 to denote the variables on the 
downstream (right) and upstream (left) sides of the shock. 
Since all the quantities on the downstream side are 
prescribed, the flow variables on the upstream side follow 
routinely from the Ranklne-Hugonlot relations. Of course 
these relations must be employed in a manner which accounts 
tor the shock velocity and curvature. 

A few preliminary definitions are needed for the equation 
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which determines the shock position as a function of the 
computational time T. Let N be the unit normal to the 

shock front. Its components in the physical plane are 

(l,3x /9y) 

N - - ^ . (9G) 

>'4 + (3x /9y)^ 

8 

Let denote the nornal velocity of a point on the shock* 

Then 


u » u N 
8 


(97) 


and Xg(Y,T) can be obtained by integrating with respect to 
T the projection of onto the X-direction* If tue 

incoming normal velocity relative to the shock is denoted by 

«rel» 

“rei ’ -“1 * ” " “s 


and the relative Mach number is 


M 


rel 


u 1 /a. • 
rel 1 


(99) 


The present numerical method presumes that ^el always 

greater than 1. 

The Ranklne-Hugonlot relations imply that 




The equation for the shock acceleration is obtained by 
differentiating Eqs* (99) and (100) and then combining the 
results to obtain 


“s.T “ ^ " 2 yM 




where 


A ■ Uj^^ • N + Uj» N^. 


( 102 ) 


The time derivatives on the right-hand-side of Eq. (101) are 
obtained from Eq. (85) using spectral approximations to the 





spatial derivatives* The shock velocity is obtained by 
integrating Bq* (101) with respect to T. 

The collocation grid in the computational plane is fixed 
and uniform. Since the shock tront moves to the right in the 
course of the calculation, the corresponding discrete grid in 
the physical plane is expanding* Thus, the effective resolu- 
tion in the x-direction continually decreases during the 
evolution. Eventually the resolution of any calculation will 
become Inadequate and the results will no longer be reliable. 
Fortunately, in many situations the important information can 
be extracted before this occurs, especially if the initial 
grid is taken to be a fine one. 

5.4. Boundary Conditions 


The correct boundary conditions at both the left and 
right boundaries depend upon the relative shock Mach number. 
If Y = 1.4 and > 2.08, then the flow behind the shock is 
supersonic. In this case both boundaries are supersonic in- 
flow boundaries and it is appropriate to prescribe all varia- 
bles there. If < 2.08, then these boundaries are subsonic 
Inflow ones. The ad' isable procedure here is to base the nu- 
merical boundary conditions on the linearized characteristics 
of the Euler equations. At the left (subsonic) boundary the 
(linearized) characteristic variable corresponding to the out- 
going characteristic direction Is 


R~ = P - - u. 
a 

similarly, 

R =■ P + - i’ 
a 


(103) 

(104) 


corresponds to the outgoing characteristic direction at the 
right (subsonic) boundary. 

A set of successful numerical boundary conditions on the 
left is obtained by first calculating preliminary values of 
all quantities at the left boundary and then incorporating the 
given values of S, v, and R’*' as 

S - S ^ 

given 


V » V . 

given 



given 


P 


Y 


It 


a 


P 


prelim 


jr 

a prelim' 


( 105 ) 
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Thus, the PDE la used to update the appropriate characLeristic 
combination of variables at the boundary* The characterlatlc 
analysis Is given in [33]* The particular numerical boundary 
condition was advocated in [34]* For the right boundary a 
similar characteristic correction procedure can be incorpora- 
ted Into the evaluation of the ?2 tena in Eq. (101)* This 
characteristic affects the shock velocity. 

At the top and bottom boundaries (which have been 
stretched to Infinity in the physical plane) zero disturbance 
boundary conditions are enforced* This is certainly justifi- 
able whenever the fluctuations decay rapidly in these direc- 
tions* However, there will be spurious reflections from the 
upper and lower boundaries If the disturbances extend that far 
out . The spurious reflections that might emanate from these 
boundaries need not pose a serious problem since the 
decreasing resolution resulting from the shock motion already 
limits the useful duration of a calculation* 

5.5. Shock Interaction Examples 

Salas , et al . [18] used the algorithm out 1 ined above to 
compute the Interaction of a shock with a single vortex, a hot 
spot and a Karman vortex street. They also gave comparisons 
with results from a similar second-order finite difference 
method. The spectral method produced virtually Identical 
results with only 1/7 as many grid points. 



Figure 6. Initial entropy contours of a 2SZ hot spot 
about to interact «ritb a Mach 1*2 shock. 

A Mach 1-2 hot spot calculation is Illustrated here in 
Figures 6 and 7* The hot spot situated In the quiescent field 
on the right in Figure 6 has the temperature distribution t 
given by 

T - k exp{-[(x-XQ)^ + (y-yQ)^]/2o^|, (1C6) 


r 


- ‘"V*! 
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where k “ 0.25, o » 1.25, = 0.5 and =* 0. The Initial 

shock position Is x » 0. Figure 7 displays the velocity 
field at time t = 0.52, after the shock wave has passed over 
the hot spot. The shock front appears a the solid line In 
both figures. 



Figure 7. Perturbed ^pstreaa velocity vectors 

after the shock-hot spot Interaction. 
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INCIDENT ANGLE 

Figure 8. Spectral (circles) and finite difference 
(squares) results for vortlclty unve 
anpllflcatlon versus Incidence angle. The 
solid line Is the linear theory prediction. 

A sample of the spectral ’•coults of Zang, et al. [30] for 
shock-turbulence Interactions Is given In Figure 8. This la a 
comparison of the computed, nonlinear amplification of Inci- 
dent Mach 3 vortlclty waves with the linear theory predlc- 
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tions* Figure 5 shows the result of a large amplitude 30° In- 
cident vortlclty wave Interacting with a Mach 1#5 shock* The 
spectral and finite difference results are comparable even 
though the spectral results were produced on a far coarser 
grid — 32 8 versus 64 ^ 32. 


Although the present spectral calculations are more effi- 
cient than the finite difference ones In terms of storage, 
they do not yet offer a clear advantage In terms of machine 
time. The culprit is the severe explicit time-step restric- 
tion for the spectral method. A robust means of surmounting 
this restriction is perhaps the most pressing need for spec- 
tral methods for evolution problems. 


5.6. The Blunt Body Problem 

The classical problem of a blunt body in a supersonic 
stream has been an Ideal test problem for numerical methods as 
it provides a relatively simple well-posed transonic problem 
with nontrivial initial and boundary conditions. Like most 
common methods the spectral method of Hussainl, et al. [19] 
obtains the steady state solution as the time asymptotic 
solution of the unsteady Euler equations which are vnrltten in 
the cylindrical polar coordinate (r,6) system. The physical 
domain of Interest consists of the known body r » r^(0), the 
unknown shock location r » rg(6,t), the axis of symmetry 
(the front stagnation streamline 0 » tt ) and the outflow 
boundary 0 « purpose of shock fitting, 

the coordinate transformation 


r - r^(0) 

8 D 

(107) 

7T - 0 
“G 

max 


Is Introduced so that the shock wave and the body are 
coordinate lines in the transformed domain. The transformed 
equations of motion, in the notation of the previous problem, 


are 


0^ + B + C Oy + R - 0. 


(108) 


where 


B » 


U 

(a^/Y)X^ 
(a^/Y)(l ,'r)XQ 


YX^ (Y/r)Xg 
U 0 

0 U 


0 

0 

0 


(109) 


0 


0 


0 


IT 



and 


with 


C 


ORiG::-^:it. is 

OF POOR QUALITY 


(a^/Y)Y 


YY^ (Y/r)Y0 
V 0 

(a‘/Y)(l/r)Y0 0 V 

0 0 0 


R- [y^. f. 0] 


0 

0 

0 

V 


( 110 ) 


( 111 ) 


+ tlX + - Xq 
t r r 9 


V - - Y„. 
r 9 


( 112 ) 


The flow field variables are expanded In double Chebyshev 
series, and the solution technique Is the same as for the 
previous problem. 

The shock boundary r « rg(9,t) (l*e., X ■ 1) Is computed 
using Ranklne-Hugonlot jump conditions and the compatibility 
equation along the Incoming characteristic from the high pres- 
sure side of the shock. At the symmetry line 9 •« tt (Y - 0) 
the 9-component of velocity v Is set equal to zero. On the 
body r * r^,(9) (l.e., X ■ 0), the normal component of 

velocity u Is zero, ^max chosen so that the outflow 
boundary Y • 1 Is supersonic, and hence no boundary 

conditions need be Imposed. 



Figure 9. Spectral solution on an 8 x 8 grid for a 

circular cylinder in a Mach 4 unlfom stream* 




I 


r 
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Olilv.l- .... F 

OF POOR QUAlin^ 


Figure 9 shows the Jfach number contours and the velocity 

vectors for a circular cylinder in a unlfo* a stream at M -4. 

00 

The results are found to be in very goa agreement with the 
tabulated values given in [36] . (The grid is so coarse that 
the contour plotter produces Jagged lines.) Figure 10 gives 
the results from the linearly-sheared stream. Even on a very 
coarse grid the spectral method captures the recirculating 
region. 
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figure 10* S|>ect>.*.'' solution on an 8 x g grid for a 

circular cylinder in a linearly-sheared streaa* 

The explicit time-step restriction is a problem here as 
well, for neither spectral solution was run to a truly 
acceptable steady-state. 


i. SamART 

Techniques are now available for obtaining viable 
spectral solutions to some compressible flow problems on grids 
far coarser than those needed for comparable finite difference 
solutions. The greatest success for shock-capturing spectral 
methods has been for potential flow. Far more sophisticated 
filtering techniques than are presently available appear 
necessary for cuccessful shock-capturing in the context of the 
Euler equations. However, when the shock is fit rather than 
captured, the Euler solutions contain no discontinuities and 
thus spectral solutions might be expected to yield exponential 
convergence. 
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